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FOREWORD 


This  report  was  prepared  by  the  University  of  Dayton  Research  Institute  under  Air  Force 
Contract  No.  F33615-95-D-5029,  Delivery  Order  No.  0004.  The  work  was  administered 
under  the  direction  of  the  Nonmetallic  Materials  Division,  Materials  and  Manufacturing 
Directorate,  Air  Force  Research  Laboratory,  Air  Force  Materiel  Command,  with  Dr.  James 
R.  McCoy  (AFRL/MLBC)  as  Project  Engineer. 

This  report  was  submitted  in  November  1999  and  covers  work  conducted  from  15 
September  1998  through  14  September  1999. 


EXECUTIVE  SUMMARY 


Three-dimensional  failure  initiation  and  progression  modeling  and  the  experimental 
work  required  to  verify  the  analysis  methods  are  reported  below.  The  stacking  sequence 
effect  on  pin  bearing  strength  was  examined  by  three-dimensional  modeling.  A  displacement 
spline  approximation  method  was  used  to  perform  the  stress  analysis.  The  laminate  and  the 
fastener  were  assumed  to  be  linear  elastic.  Processing  residual  stresses  were  included  in  the 
analysis.  The  contact  region  and  stress  were  unknown  a  priori  and  resulted  from  the 
solution.  A  nonuniform  contact  region  through  the  thickness  was  allowed.  Undamaged 
thermoelastic  moduli  were  employed  in  the  model.  The  stacking  sequence  effect  on  bearing 
strength  was  investigated  by  comparing  stress  distributions  in  two  quasi-isotropic  laminates 
of  [02/902/452/-452]s  and  [02/452/902/-452]s  stacking  sequence.  Stress  magnitudes  were 
displayed  relative  to  ply  strength  properties.  It  was  observed  that  all  stress  components, 
including  the  in-plane  stresses,  were  affected  by  the  stacking  sequence.  In  the  vicinity  of  the 
bearing  plane,  very  high  transverse  shear  stresses  were  found.  Using  a  point  stress  failure 
criterion,  a  prediction  of  failure  onset  was  performed.  The  predicted  failure  initiation  load 
for  the  [02/902/452/-452]s  stacking  sequence  was  approximately  40  percent  higher  than  the 
predicted  failure  initiation  load  for  the  [02/452/902/-452]s  stacking  sequence. 

Damage  initiation  and  progression  in  the  unidirectional  and  quasi-isotropic  laminates 
was  modeled  based  on  the  property  degradation  methodology.  The  Reissner-Hellinger 
variational  principle  was  applied  to  rigorously  represent  the  degraded  material  properties  in 
order  to  describe  discrete  damage  at  the  hole  edge.  A  recursive  algorithm  based  on  the 
damage  increment  concept  was  developed  to  perform  the  damage  evolution  modeling.  Both 
mechanical  loading  and  residual  stresses  were  taken  into  account.  The  maximum  stress 
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failure  criterion  was  implemented  and  used  to  select  the  elements  where  the  properties  were 
degraded.  Conventional  finite  element  formulation  based  on  linear  displacement 
approximation  was  used.  The  analysis  uncovered  severe  mesh-type  dependence  of  the 
damage  evolution  modeling  based  on  conventional  finite  element  approximation.  The  results 
showed  that  only  a  specific  class  of  finite  element  meshes  allows  one  to  describe  the  stress 
relaxation  in  the  fiber  direction  in  unidirectional  laminates  with  holes  occurring  due  to  matrix 
splitting  under  uniaxial  tension.  Further  research  will  be  conducted  to  employ 
unconventional  methods  such  as  spline  mesh  overlays  and  spline  basis  enrichment  approach 
to  discrete  damage  modeling  in  laminates  with  holes. 

An  experimental  investigation  has  been  conducted  on  the  initiation  and  growth  of 
damage  in  close  proximity  to  open  hole  for  the  [±30/902] sand  [02/902]s  laminates  under 
incremental  tension  loading.  For  both  laminates,  the  analytical  predictions  of  strain 
calculated  by  spline  variational  theory  and  by  elasticity  showed  good  agreement  with 
experimentally  measured  results.  The  development  of  damage  observed  by  x-radiography 
was  correlated  with  the  stress-strain  behavior.  As  damage  developed  for  both  laminates, 
stress  concentration  decreased,  and  the  net  stress  at  failure  was  nearly  identical  to  the  static 
strength  of  the  respective  unnotched  coupon.  However,  for  the  first-ply  failure,  the  net  stress 
onset  is  related  through  a  reduced  stress  concentration  factor.  This  experimental  work 
continues  on  the  [±45/90/0]  s  and  [45/0/-45/90]  s  laminates. 
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1.  THREE-DIMENSIONAL  STRESS  ANALYSIS  AND  FAILURE  PREDICTION 

IN  FILLED-HOLE  LAMINATES 


Bolted  joining  remains  a  primary  joining  method  of  load  carrying  composite 
structural  elements  in  the  modem  aerospace  industry.  Significant  attention  in  the  literature 
has  been  given  to  practically  all  aspects  of  composite  bolted  joining,  both  in  the  analysis  and 
experimental  fields  as  described  in  a  recent  survey  [1] ,  The  need  for  a  composite  bolted  joint 
design  tool  led  to  development  of  several  composite  bolted  joint  failure  prediction  programs, 
reviewed  by  Snyder  et  al.  [2] ,  and  utilized  in  the  industry.  The  common  mechanical 
foundation  for  all  of  these  tools  reviewed  is  lamination  theory  combined  with  average  or 
point  stress  failure  criterion  to  predict  final  strength.  A  different  failure  prediction 
methodology,  based  on  fracture  mechanics,  was  developed  in  [3,4]  and  is  based  on  pseudo¬ 
crack  propagation  in  the  radial  direction  from  the  hole  edge.  Energy  release  rates  for  these 
cracks  are  obtained  by  utilizing  the  rule  of  mixtures  and  the  number  of  plies  of  each 
orientation  -  0,  ±45,  and  90  -  in  the  laminate.  In  all  of  the  above  methods,  the  fastener-plate 
contact  stress  is  assumed  to  vary  as  a  cosine  function  over  half  the  hole  circumference  with 
an  amplitude  of  -4/71  times  the  nominal  bearing  stress.  The  advantage  of  these  largely 
empirical  approaches  is  that  a  broad  spectrum  of  problems  can  be  solved  without  significant 
computer  resources.  The  disadvantage  is  the  need  for  an  extensive  test  database  that  is 
required  for  each  new  material  system  utilized.  To  develop  new  tools  for  composite  bolted 
joint  design  and  strength  prediction  based  on  measured  ply  properties,  an  in-depth 
understanding  of  the  bearing  failure  mechanisms  is  important.  Recent  experimental  works 
[5-7]  dealing  with  graphite-epoxy  laminates  provide  detailed  investigations  of  the  bearing 
failure  process  and  show  similar  results.  The  specimens  subjected  to  bearing  loading  were 
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x-rayed  and  sectioned  at  different  load  levels  lower  or  equal  to  the  failure  load.  Microscopic 
studies  of  the  bearing  plane  sections  showed  that  bearing  failure  consists  of  a  damage 
accumulation  process  initiated  in  the  interior  of  the  laminate  at  the  hole  edge.  It  consists  of 
45°  transverse  shear  cracks  propagating  to  the  outer  surfaces  from  the  hole  edge.  According 
to  [6,7]  in  the  case  of  a  pin-loaded  hole,  these  cracks  reach  the  outer  surfaces  at 
approximately  a  quarter  of  the  laminate  thickness  away  from  the  hole  edge,  resulting  in  final 
failure.  In  all  cases  massive  delamination  was  observed  on  the  interior  interfaces  prior  to 
failure.  In  the  presence  of  clamping  forces  [8],  the  damage  mechanism  does  not  seem  to 
change.  The  failure  load  significantly  increases  because  the  final  failure  does  not  occur  until 
the  damage  accumulation  region  reaches  the  outer  washer  diameter,  provided  that  the  washer 
can  carry  the  increasing  plate  expansion  force.  As  concluded  in  [5,6] ,  the  bearing  failure  is  a 
substantially  three-dimensional  process  involving  delamination  and  transverse  shear; 
although  in  the  sequel  of  works  [7,8] ,  the  authors  proposed  a  two-dimensional  bearing 
strength  prediction  model,  which  was  shown  to  provide  good  correlation  with 
experimentally-measured  bearing  strength.  This  model  is  based  on  nonlinear  finite  element 
simulation  with  gradual  property  degradation  algorithms. 

Lamination  theory,  however,  cannot  account  for  stacking  sequence  effects  on  bearing 
strength.  Hamada  et  al.  [5]  tested  four  symmetric  laminates  stacked  with  sequences  of  O2, 
452,  -452,  and  902  ply  groups.  The  material  used  was  T300/2500.  The  bearing  strength 
results  obtained  in  [5]  are  summarized  in  Table  1,  and  the  specimen  geometry  used 
throughout  this  research  is  shown  in  Figure  1.  An  approximately  13  percent  difference  was 
observed  for  specimens  C  and  N,  corresponding  to  stacking  sequences  [02/902/452/-452]s  and 
[02/452/902/-452]s,  respectively.  These  specimens  differ  by  an  apparently  minor  alteration  in 
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Table  1 

Experimental  Bearing  Failure  Stresses  Obtained  from  [1] 


Specimen  Identification  and 

Stacking  Sequence 

N  -  [02/+452/902/-452]s 

318 

A  -  [02/+452/-452/902]s 

333 

B  -  [-452/02/+452/902]s 

338 

C  -  [02/902/+452/-452]s 

360 

Figure  1.  Specimen  Geometry.  L=72  mm,  A=30  mm,  D=6  mm,  Xc=  12  mm,  Yc=15  mm, 
and  ply  thickness  ~0.25  mm. 
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stacking  sequence  only.  The  only  difference  is  in  exchanging  the  90°  and  45°  plies  while 
retaining  the  position  of  the  0°  plies  on  the  outside  and  -45°  plies  at  the  midsurface.  The 
purpose  of  the  present  work  is  to  perform  three-dimensional  stress  analysis  to  explain  the 
effect  that  such  a  small  change  of  the  stacking  sequence  can  have  on  laminate  pin  bearing 
strength.  Damage  related  stress  redistribution  was  not  considered  in  this  research. 

A  previously  developed  fully  three-dimensional  stress  analysis  procedure  [9]  was 
used  to  examine  the  stress  distributions  occurring  under  bearing  loading.  The  procedure  is 
based  on  a  cubic  spline  approximation  of  displacements  along  with  a  curvilinear 
transformation  technique  also  based  on  the  spline  approximation.  Due  to  the  highly 
nonuniform  stress-strain  field  through  the  thickness  of  the  laminate,  the  contact  zone 
becomes  nonuniform  through  the  thickness.  An  automated  algorithm  was  developed  for 
contact  zone  definition.  Verification  of  the  results  of  the  interlaminar  stress  prediction  was 
accomplished  by  comparison  with  the  asymptotic  solution. 

1.1  Problem  Definition 

Consider  a  rectangular  orthotropic  plate  containing  a  circular  hole  having  a  diameter 
D,  as  shown  in  Figure  1.  The  plate  consists  of  N  plies  of  total  thickness  FI  in  the  z-direction, 
has  a  length  L  in  the  x-direction  and  width  A  in  the  y-direction.  A  circular  isotropic  elastic 
fastener  of  diameter  d  is  situated  at  the  center  of  the  hole.  The  length  of  the  fastener  in  the 
z-direction  is  equal  to  the  plate  thickness  H.  We  will  consider  the  bearing  loading  case  when 
the  load  is  applied  via  displacement  boundary  conditions  at  x=L,  while  all  other  boundaries 
were  kept  free: 


ux{L,y,  z)  =  uL,  uv(L,y,z)  =  u.(L,y,z )  =  0, 
<7n.  (0,  y,  z)  =  <J  r (0 ,y,z)  =  crr.(0,  y,  z)  =  0 
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A  cylindrical  coordinate  system  is  defined  originating  from  the  center  of  the  hole: 

x  =  rcos©  +  xc,  y=  /-sin  ©+  yc,  z  =  z  (2) 


where  Xc  and  yc  are  the  coordinates  of  the  center  of  the  hole.  The  fastener  displacements  will 
be  denoted  as  vr,  ve,  and  vz.  The  in-plane  motion  of  the  fastener  is  fixed  at  two  points  located 
at  the  center  of  the  top  and  bottom  fastener  surfaces  z=0  and  z=H,  and  the  z  displacement  is 
fixed  at  the  mid-surface.  In  cylindrical  coordinates  the  boundary  conditions  are  as  follows: 


vr(o,e,o)=vfl(o,e,o)=o, 

v,(0,©,/7)  =  ve(0,©,AO  =  0,  (3) 

v C(O,0,///  2)~  0 

For  small  deformations,  the  nonpenetration  condition  can  be  simplified  as 

ur{DI2,d,z)-vr{dl2,d,z)  +  M?>U  ,  (4) 

where  AR—D/2-ci/2  is  the  clearance  between  the  hole  and  the  fastener,  and  u,  is  the  radial 
displacement  of  the  composite  plate. 

The  minimum  potential  energy  principle  along  with  the  Lagrangian  multiplier  method 
were  used  to  solve  the  problem.  Frictionless  contact  is  considered.  The  first  variation  of  the 
following  functional  is  required  to  vanish. 


n  +n  + 

P  b 


JI*«  z)(ur(DI2Az)-vr(d/2A  z)+M2)dS 


(5) 


V  ate,-)  ) 

The  potential  energy  of  the  laminate  is  denoted  as  ITp,  and  the  potential  energy  of  the  elastic 
fastener  as  nb.  £2(0, z)  is  the  contact  zone  on  which  the  relationship  (4)  becomes  an 
equality.  The  size  and  location  of  this  zone  is  unknown  initially.  Functional  (5)  suggests  a 
simple  interpretation  for  the  Lagrangian  multiplier  as  the  value  of  the  radial  stress  at  the 
contact  surface,  i.e.,  (D  /  2,0  ,z)  =  oA  (d  /  2,0  ,z)  =  A(0,z);(0,z)e  0(0, z). 
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In  the  following  sections,  the  spline  approximation  of  the  functions  included  in  (5) 
will  be  constmcted. 

1.2  Spline  Approximation 

Curvilinear  coordinates  p  and  (j)  were  introduced  to  map  the  plate  with  a  cutout  into  a 
region  0  <  p<  1  and  0  <  (j)  <  2tt  .  Coordinate  lines  of  this  transformation  are  shown  in  Figure  1. 
The  transformation  is  built  according  to  [10]  so  that  the  coordinate  line  p=0  corresponds  to 
the  hole  edge  and  p=l  to  the  outer  rectangular  contour  of  the  plate.  The  radial  lines 
correspond  to  (^constant.  Sets  of  B-type  cubic  basis  spline  functions 


along  each  coordinate  were  built  upon  subdivisions  0  =  p0  <p,  <  ...<  pm  =  1, 

0  =  (p0  <  0,  < ...  <  <4  =  2 n,  and  C'1'1  =  zn  <  q  < ...  <  z„  =  z01'  so  that  the  s-th  ply  occupies  a 
region  ]s~l)  <  z<  z(s),  and  n,  is  the  number  of  sublayers  in  each  ply.  The  subdivision  of  the  p 
coordinate  is  essentially  nonuniform.  The  interval  size  increases  in  geometric  progression 
beginning  at  the  hole  edge.  The  region  0  <p  <ph  in  which  the  curvilinear  transformation  is 
quasi-cylindrical  is  subdivided  into  rr^  intervals,  so  that  ph  =  pniQ .  The  number  of  intervals 

of  subdivision  m,  k,  and  ry  in  each  direction  along  with  the  mesh  nonuniformity 
characteristics,  such  as  rpj  and  the  consecutive  interval  ratio,  determine  the  accuracy  of  the 
solution  and  the  size  of  the  problem. 

The  detailed  solution  procedure,  including  determination  of  the  contact  region,  has 
been  described  in  [9] . 


8 


1.3  Failure  Criteria 


Bearing  failure  can  be  manifested  by  a  combined  effect  of  all  individual  failure 
modes:  tensile  fiber  breakage,  micro-buckling,  delamination,  matrix  cracking,  and 
compression/shear  cracking.  As  a  noncatastrophic  failure  mode  [6] ,  it  can  only  be 
adequately  described  by  taking  into  account  the  stress  redistribution  resulting  from  damage 
initiation  and  growth.  As  a  result  of  the  three-dimensional  nature  of  the  stress  field  in  the 
bearing  contact  region,  the  modeling  of  bearing  failure  represents  a  challenging  task.  It  has 
been  accomplished  so  far  in  simplified  two-dimensional  models  [7,8]  that  neglect 
delamination  and  stacking  sequence  effects.  In  the  present  work,  a  three-dimensional  stress 
analysis  was  performed  instead  to  evaluate  the  magnitude  of  the  interlaminar  stress 
components  and  the  effect  of  stacking  sequence  on  bearing  strength.  No  effect  of  stress 
redistribution  due  to  damage  was  considered. 

To  put  the  stress  values  in  perspective,  they  were  compared  to  the  respective  ultimate 
strength.  Six  dimensionless  functions  s\ ,  6, were  introduced  such  that 


^rr,o n  >o 

Jl  • -p 

ii  <  o 

.  Ac- 


Y~^22  >  0  y-^33  >  0 

~L’<J22  <  0  <0 


<*23 

t’  HZ 

<7,3 

c  = 

^12 

5 

^5 

5 

^6 

5 

(7) 


(8) 


(9) 


The  stress  components  Cij,  i,j=  1,2,3  are  in  ply  material  coordinates,  and  the  xi-axis  coincides 
with  the  fiber  direction.  Parameters  XT  and  Xc  are  the  tension  and  compression  strength  of 
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the  unidirectional  laminate  in  the  fiber  direction,  Yt  and  Yc  are  the  tension  and  compression 
strength  in  the  transverse  direction,  and  S  is  the  in-plane  shear  strength.  The  in-plane  shear 
strength  is  also  used  for  scaling  the  interlaminar  shear  stress  components. 

Bearing  loading  generates  a  three-dimensional  state  of  stress  in  the  vicinity  of  the 
hole  edge.  Comparison  of  the  stress  components  to  the  respective  ultimate  material  direction 
strengths  neglects  mode  interaction.  In  the  present  work,  the  use  of  interactive  failure 
criteria,  such  as  Tsai-Wu  [11]  or  Hashin  [12] ,  has  been  avoided  because  of  a  lack  of 
experimental  evidence  supporting  the  reliability  of  strength  prediction  by  these  criteria  in  a 
three-dimensional  state  of  stress  as  well  as  under  stress  gradients.  All  stress  components  in 
the  present  paper  will  be  examined  separately  by  using  Equations  (7)-(9). 

Three-dimensional  ply  level  analysis  performed  in  the  present  paper  results  in  stress 
singularities  at  the  ply  interface  and  the  hole  edge  intersections.  In  order  to  predict  the  onset 
of  failure,  a  point  stress  criterion  was  applied.  Failure  initiation  was  predicted  when  one  of 
the  functions  s;,  i=l,...,6,  has  an  absolute  value  greater  than  one  at  a  distance  of  a  ply-group 
(two  plies)  thickness  away  from  the  hole  edge.  The  choice  of  a  one  ply-group  thickness 
characteristic  distance  was  influenced  by  Kim  and  Soni  [12]  who  used  a  one-ply  thickness 
characteristic  distance  with  an  average  stress  failure  criteria  for  prediction  of  delamination 
onset. 

1.4  Numerical  Results  and  Discussion 

A  schematic  of  the  bearing  specimen  is  shown  in  Figure  1.  To  perform  the 
dimensional-dimensional  stress  analysis  and  to  apply  the  failure  criteria,  a  full  set  of 
thermoelastic  and  strength  characteristics  is  needed.  Because  of  a  lack  of  sufficient  data  on 
the  material  tested  in  [6] ,  a  similar  T300/5250  material  was  used  for  the  analysis.  The 
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stiffness  and  strength  parameters  for  T300/5250  are  given  in  Table  2.  The  estimated 
temperature  difference  between  the  room  temperature  and  the  cure  temperature,  required  for 
residual  stress  calculation,  was  -100°C.  Bearing  load  was  applied  through  the  displacement 
boundary  conditions  (1).  The  applied  bearing  stress  was  calculated  as 

j  H  A 

1  y>  *)  dydz  ( 1  o) 

ADi  t 


Table  2 

Material  Properties  Used  in  Modeling 


Elastic  Constants 

Failure  Parameter 

Ei  =  181.0  GPa 

E2  =  E3  =  10.3  GPa 

V  ]  2  =  V  13  =  0.28 

V23  =  0.40 

Gn  =  Gi3  =  7.17  GPa 

G23  =  3.68  GPa 
a,  -  0.02x1 0“6/°C 
a2  =22.5x1 0'6/°C 

XT=  1500  MPa 

Xc=  1500  MPa 

Yt  =  40  MPa 

Yc  =  246  MPa 

S  =  68  MPa 

The  subdivisions  in  the  circumferential  and  radial  directions  in  the  plate  were  k=36 
and  m=18,  with  mo=12  intervals  in  the  near  hole  region  and  the  size  ratio  q=l  .2.  The 
subdivisions  in  the  radial  and  circumferential  directions  in  the  fastener  were  kb=36,  mb=8, 
and  qb=1.2.  Two  sublayers  through  the  thickness  of  each  ply  were  used.  Only  the  upper  half 
of  the  laminate  was  modeled  because  of  symmetry. 

First,  the  effect  of  stacking  sequence  on  stress  components  was  examined.  The  load 
level  corresponding  to  failure  of  the  weaker  of  the  two  specimens  was  chosen,  so  that 
Ob=3  1 8MPa.  The  influence  of  stacking  sequence  on  the  in-plane  stress  components  is 
illustrated  in  Figure  2.  The  functions  .s’] ,  S2,  and  .sy„  introduced  in  Equations  (7)-(9),  are 
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S2 - specimen  N:  [02/+452/902/-452]s 


function  value  scale 
-1,0  0.0  1.0 


Figure  2.  Stacking  Sequence  Effect  on  the  In-Plane  Stress  Components  at  the 
Bearing  Plane  for  Both  Specimens. 


12 


shown  on  the  bearing  plane.  This  location  was  chosen  due  to  available  experimental 
evidence  showing  the  bearing  plane  to  be  a  critical  damage  area  in  the  specimens.  Both 
specimens  had  holes  with  a  relatively  small  diameter  (compared  to  specimen  dimensions) 
and  exhibited  bearing-type  failure.  Five  gray-scale  color  codes  are  used  to  display  the 
results.  The  white  regions  correspond  to  .s-,<- 1  (i.e.,  compression  failures),  and  the  darkest 
regions  correspond  to  S{>1  (i.e.,  tensile  failures)  for  i=l,2,3.  In  the  case  of  the  shear  stresses, 
the  functions  y,  1=4,5, 6,  can  only  have  positive  values.  The  largest  effect  of  the  stacking 
sequence  change  on  in-plane  stress  components  can  be  seen  for  the  transverse  normal 
stresses.  Figure  2b  shows  that  the  .s’2>l  region  (tensile  matrix  cracking)  in  the  0°  ply  is  twice 
as  large  in  specimen  N  compared  with  specimen  C.  This  can  be  explained  by  the  strong 
reinforcing  effect  of  the  90°  ply  in  specimen  C  compared  with  the  weaker  reinforcing  effect 
of  the  45°  ply  in  specimen  N.  The  fiber  direction  normal  and  in-plane  shear  stresses,  Figures 
2a  and  2c,  do  not  show  a  significant  stacking  sequence  effect.  It  should  be  mentioned  that 
the  function  values  should  be  compared  in  the  same  plies  in  the  two  specimens.  For 
example,  the  difference  in  appearance  of  the  shear  stress  function,  %  in  Figure  2c  is 
deceptive.  The  distribution  of  sc,  in  the  45°  and  90°  plies  in  both  specimens  is  quite  similar. 

Interlaminar  stresses  developed  at  the  bearing  plane  were  also  examined.  Of  the  three 
interlaminar  components,  the  transverse  shear  stress,  Ob,  was  found  to  extend  over  the 
largest  region  where  the  ultimate  value  was  exceeded.  Figure  3c  shows  the  contour  plots  of 
the  function  S5  [Equation  (9)].  For  both  specimens,  Ob  exceeded  its  ultimate  value  in  the 
vicinity  of  the  hole  edge.  In  the  case  of  the  stronger  specimen,  C,  the  area  where  is 
much  smaller  and  consists  of  two  isolated  subregions  at  the  90/45  and  45/-45  interfaces.  The 
other  interlaminar  shear  stress  component,  O23,  reflected  in  Figure  3b  by  function  sq,  exhibits 
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a) 


rzu 


S3 - specimen  C:  [02/902/+452/-452]s 


b) 


c) 


S5 - specimen  C:  [02 /9 02/+452 /-4 5 2 ] s 


S5 - specimen  N:  [02/+452/902/-452]s 


Figure  3. 


function  value  scale 
-1,0  0.0  1.0 


Stacking  Sequence  Effect  on  the  Interlaminar  Stress  Components  at  the 
Bearing  Plane  for  Both  Specimens. 
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lower  values.  The  function  .sy.  related  to  the  transverse  normal  stress,  is  shown  in  Figure  3a. 
Compared  to  s$  it  also  displays  lower  values;  however,  a  significant  difference  in  the 
distributions  between  the  two  laminates  can  be  seen.  The  weaker  laminate,  N,  displays 
significant  tensile  stresses  throughout  the  laminate  thickness,  whereas  laminate  C  has 
significant  transverse  normal  tensile  stress  only  in  a  small  region  in  the  90°  ply  and  at  the 
90/45  interface.  This  is  in  qualitative  agreement  with  the  cross-sectional  micrographic 
studies  of  the  failed  specimens  in  Reference  [6] .  The  cross  section  of  specimen  N  shows 
significant  expansion  in  the  z-direction  with  open  delaminations  in  addition  to  transverse 
shear-induced  45°  matrix  cracking  and  fiber  buckling  bands.  In  the  case  of  specimen  C,  the 
failed  cross  section  is  significantly  less  expanded  in  the  transverse  direction,  with 
delamination  clearly  visible,  but  with  much  less  opening.  It  must  be  noted  that  the 
correlation  outlined  above  is  qualitative,  because  the  stress  state  obtained  under  the 
assumption  of  a  virgin  specimen  is  not  valid  at  the  advanced  stages  of  damage. 

The  x-rays  of  the  two  failed  specimens,  obtained  from  [6] ,  show  different  patterns  of 
0°  ply  splits.  The  failure  load  was  Gb  =318  MPa  and  Gb  =360  MPa  for  specimen  N  and 
specimen  C,  respectively.  The  functions  si  and  S(,  were  analyzed  at  these  load  levels  in  the  0° 
ply  in  both  specimens  near  the  interface  with  the  underlying  ply.  The  observed  matrix 
cracking  can  be  produced  both  by  the  in-plane  transverse  tensile  stress  component  and  the 
in-plane  shear  stress  component.  Figure  4  displays  superimposed  gray-scale  plots  of  the 
functions  .V2  and  ,sy,.  Only  two  .v,  value  intervals,  ,s-i<  1  and  .v,>  l ,  are  shown  for  each  function. 
The  white  area  corresponds  to  ,v,<  1 .  The  light  gray  corresponds  to  .sy,>  1 ,  the  dark  gray  to 
52>1,  and  the  medium  gray  designates  the  area  where  both  $(>  1 .  Figure  4  also  displays  the 
schematics  of  the  matrix  cracks  in  the  0°  plies.  They  are  distributed  in  a  relatively 
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Specimen  C:  [02/902/+452/-452]s 


Specimen  N:  [02/+452/902/-452]s 


Sg>l  S2>1  &  Sg>l  So>l 


Figure  4.  Comparison  of  s2  and  s6  for  0°  Plies  in  Both  Specimens  at  their 

Respective  Failure  Loads.  The  straight  solid  lines  represent  the  most 
clearly  distinguishable  0°  ply  cracks  obtained  from  x-rays  published  in 
Reference  [1]. 


symmetrical  manner  in  laminate  C  and  in  a  distinctly  non  symmetrical  pattern  in  laminate  N. 
These  crack  distributions  agree  well  with  distributions  of  the  combined  fields  of  52  >1  and 
i’6>  1  in  each  laminate  displayed  in  Figure  4.  It  is  worth  noting  that  the  application  of 
Hashin’s  [12]  matrix  cracking  criteria,  which  combines  the  shear,  transverse  normal,  and 
interlaminar  normal  stress  components,  also  provides  good  agreement  with  the  0°ply 
splitting  patterns  in  the  two  specimens. 

Failure  initiation  prediction  was  performed  by  examining  criteria  (7)-(9)  at  one  ply- 
group  thickness  away  from  the  hole  edge.  The  two  specimens  appeared  to  have  the  same 
initial  damage  mode  -  tensile  matrix  cracking  in  the  0°and  90°  plies.  Figure  5  shows  the  52 
distribution  in  these  plies  at  the  load  levels  when  the  .s'2>l  region  reached  the  dashed  contour, 
which  is  a  double  ply  thickness  (one  ply-group)  away  from  the  hole  edge.  This  occurred  in 
the  0°  ply  at  80  MPa  and  in  the  90°  ply  at  96  MPa  in  specimen  N  and  in  the  0°  ply  at  128 
MPa  and  in  the  90°  ply  at  1 12  MPa  in  specimen  C.  In  the  0°  ply  this  region  is  in  the  vicinity 
of  the  bearing  plane,  and  in  the  90°  ply  there  are  two  regions  at  opposite  sides  of  the  hole. 
Interestingly,  the  cracking  of  0°  and  90°  plies  is  predicted  in  reverse  order  for  the  two 
laminates.  The  tensile  transverse  stress  component  is  the  dominating  stress  component  at 
these  locations. 

1.5  Conclusions 

1.  Three-dimensional  stress  analysis  in  quasi-isotropic  laminates  under  bearing  loading  was 
performed  to  explain  the  effect  of  stacking  sequence  on  the  pin  bearing  strength  observed 
in  [6]  by  calculating  accurate  stress  fields  in  undamaged  plies. 

2.  A  significant  effect  of  stacking  sequence  was  observed  on  transverse  normal  in-plane 
stress  component  and  predicted  failure  initiation  loads.  The  specimen  with  the  lower 
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0°  Ply  at  128  MPa  bearing  stress 
on  the  top  face  of  the  laminate 


near  the  90/+45  interface 


Specimen  C:  [(>2/902/+452/-452]s 


0°  Ply  at  80  MPa  bearing  stress  90°  ply  at  96  MPa  bearing  stress 

on  the  top  face  of  the  laminate  near  the  90/-45  interface 

Specimen  N:  [Q2/+452/902/-452]s 


function  value  scale 


-LG  0,0  1,0 


Figure  5.  Initial  Damage  Prediction  for  Both  Specimens  (Function  s2).  The  dashed 
lines  represent  a  radius  equal  to  the  hole  radius  plus  1  ply-group 
thickness. 
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bearing  strength  exhibited  lower  damage  initiation  loads  and  a  larger  critical  transverse 
normal  tensile  stress  ratio  area  in  plies  with  the  same  orientation.  The  0°  ply  cracking 
patterns  observed  in  experiments  agree  with  superimposed  in-plane  shear  and  transverse 
normal  stress  maps. 

3.  Transverse  shear  and  normal  stresses  were  analyzed  at  the  bearing  plane  of  the 
specimens.  A  larger  area  of  critical  transverse  shear  stress  ratio  was  found  in  the 
specimen  with  lower  bearing  strength.  The  transverse  normal  stress  at  the  bearing  plane 
of  the  weaker  specimen  had  large  areas  of  significant  tensile  values  correlating  with 
massive  thickness  expansion  and  delamination  under  bearing  loading.  The  stronger 
specimen  had  smaller  areas  of  significant  tensile  transverse  normal  stresses  through  the 
thickness  of  the  bearing  plane,  which  is  consistent  with  the  smaller  thickness  expansion 
and  mostly  compression  shear  failure  pattern  observed. 
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2.  PROPERTY  DEGRADATION  BASED  PROGRESSIVE  DAMAGE  MODELING 
IN  OPEN-HOLE  COMPOSITE  LAMINATES 


To  provide  accurate  strength  prediction  in  notched  composites,  it  is  important  to  take 
into  account  the  stress  redistribution  due  to  damage  progression  prior  to  failure.  Finite- 
element  formulation  and  property  degradation  rules  are  used  through  [14-26]  to  address  this 
issue.  Two-dimensional  property  degradation-based  models  were  developed  for  strength 
prediction  in  open-hole  [14-17],  pin-loaded  [18]  and  bolted  with  clamp-up  [19-21] 
composite  laminates.  The  effect  of  clamp-up  in  [19-21] ,  which  is  essentially  a  three- 
dimensional  effect,  was  taken  into  account  by  introducing  a  washer-plate  friction  correction 
to  the  bearing  load.  This  allowed  estimation  of  the  clamp-up  effect  within  two-dimensional 
formulations,  based  on  lamination  theory.  Papers  [14-16,18]  used  essentially  similar 
property  degradation  models,  consisting  of  Hashin’s  [12]  failure  criteria  and  selected  elastic 
moduli  discounts  within  the  element.  The  thermal  residual  stresses  were  neglected.  Both 
linear  and  nonlinear  geometric  and  constitutive  modeling  were  considered  in  [18].  Five  local 
failure  options  were  implemented:  (i)  matrix  tensile  failure  resulting  in  zero  transverse 
modulus,  (ii)  matrix  compressive  failure  also  resulting  in  zero  transverse  modulus,  (iii)  fiber- 
matrix  shearing  failure  resulting  in  zero  in-plane  shear  modulus,  (iv)  fiber  tensile  failure 
resulting  in  complete  stiffness  loss  of  the  element,  and  (v)  fiber  compression  (bearing) 
failure.  Different  post-failure  property  degradation  rales  were  used  for  this  (v)  failure  mode 
in  the  linear  and  nonlinear  analysis  models.  In  the  linear  analysis  the  effect  of  fiber 
compression  failure  was  the  same  as  fiber  tensile  failure  and  resulted  in  complete  loss  of 
element  stiffness.  On  the  contrary,  in  the  nonlinear  analysis  model  the  failure  mode  (v) 
resulted  in  stiffening  the  element  by  increasing  all  elastic  moduli  10  times.  The  nonlinear 


20 


model  included  geometric  nonlinearity  and  cubic  in-plane  shear  stress-strain  relationship 
with  nonlinearity  coefficient,  which  is  obtained  by  best  fit  of  the  [±45]s  coupon  stress-strain 
curves.  It  was  concluded  that  the  final  failure  load  can  be  predicted  usually  to  within  20 
percent  of  the  experimental  results.  Inclusion  of  the  nonlinear  effects  was  found  to  improve 
the  predictions.  It  is  of  interest,  however,  as  to  how  much  of  that  improvement  is  due  to 
changing  the  property  degradation  rule  for  the  fiber  compression  failure  mode  in  the  critical 
bearing  failure  region  from  loss  of  stiffness  to  “incompressibility.” 

A  discrete  internal  state  variable  approach  to  property  degradation  was  taken  in  [1 7], 
Two  types  of  damage,  matrix  failure  and  fiber  breakage,  were  considered.  The  stiffness 
degradation  in  the  lamina  did  not  depend  on  the  cause  of  the  damage.  A  quadratic,  non¬ 
interactive  fiber  failure  criterion  was  applied.  Matrix  cracking  was  considered  to  occur  in  the 
element  if  the  Tsai-Wu  [11]  criterion  was  satisfied,  and  the  fiber  failure  criterion  was  not 
met.  The  degraded  stiffness  moduli  were  represented  as 

F  =  D  F°  F  =  D  F°  G  =  D  G° 

^11  1 1  ?  J-^22  l^2l~‘22'>  w12  Iy6u12’ 

Empirical  knockdown  factors  Dj,  Dj  and  Dr,  were  used  to  degrade  the  stiffness  and,  as 
expected,  the  values  of  these  parameters  had  a  significant  effect  on  final  failure  prediction. 
The  transverse  cracking  related  factors  were  constant  for  all  cases  Dr  =Dj  =0.2,  and  the  fiber 
breakage  related  factor  Dj  varied  from  0.07  for  the  AS4/3502  system  to  0.001  for  the 
T300/1034-C  system.  Thermal  and  moisture  effects  were  incorporated.  Interestingly,  the 
residual  stresses  in  all  cases  resulted  in  higher  predicted  strength  values.  In  the  damage 
accumulation  studies  reported,  the  model  tended  to  under-predict  the  fiber  breakage  initiation 
loads  and  over-predict  the  extent  of  fiber  failure  regions.  This  likely  resulted  from 
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difficulties  in  accurately  describing  the  local  fiber  stress  relaxation,  due  to  transverse  splitting 
and  smoothing  the  drastic  effect  of  stress  redistribution  due  to  fiber  failure,  by  considering 
only  0.001  to  0.07  times  knockdown  factors,  respectively,  for  this  failure  mode.  These 
effects  can  be  accurately  predicted  only  by  using  3-D  modeling  methods.  With  a  proper 
choice  of  the  property  degradation  factors,  the  strength  of  laminates  reported  in  the  paper  was 
predicted  to  within  14  percent. 

An  internal  state,  variable-matrix  crack  density  was  introduced  in  [21]  to  modify  the 
nonlinear  property  degradation  models  of  [14-16,18].  Exponential,  property  degradation 
rules  resulted  in  more  stable  numerical  mechanisms.  It  was  concluded  in  [22]  that  the 
developed  method  is  capable  of  predicting  the  strength  and  failure  mode  of  composite  joints 
under  combined  bearing  and  axial  bypass  loads.  It  was  also  noted  that  to  fully  account  for 
the  factors  affecting  the  strength  of  composite  bolted  joints,  the  analysis  needs  to  be  extended 
to  three  dimensions. 

Recently  several  three-dimensional  progressive  damage  models  were  developed  for 
strength  prediction  in  open-hole  [23,24]  and  bolt-hole  laminates  [25,26],  Three-dimensional 
property  degradation  modeling  described  in  [23]  was  facilitated  in  an  ABAQUS  environment 
and  consisted  of  modified  Hill  failure  criteria  described  in  [18]  and  associated  property 
degradation  rules.  The  following  damage  modes  were  incorporated:  fiber  failure, 
corresponding  to  property  degradation  rule  En=Gi2=0,  Vi3=vi2=0,  matrix  in-plane  failure, 
corresponding  to  E22=Gi2=0,  V3i=V32=0,  and  matrix  failure  due  to  out-of-plane  stress 
components,  modeled  as  E33=G3i=  G32=0,  V2]=V23=0.  To  model  the  delamination  within  the 
property  degradation  framework,  the  authors  introduced  a  resin  reach  layer  with  the 
properties  of  cured  matrix  between  each  ply.  Due  to  computational  constraints  this  layer  had 
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the  same  thickness  as  the  unidirectional  plies.  Several  examples  were  considered.  Axial 
splits  in  a  0°  ply  emanating  from  a  crack-type  notch  under  axial  tension  were  modeled.  Good 
agreement  of  the  experimentally  observed  split  length  as  a  function  of  load  was  shown. 
Symmetric  cross-ply  laminates  with  open  holes  under  uniaxial  tension  were  also  considered. 
It  was  stated  that  the  strength  of  the  laminate  was  predicted  with  five  percent  accuracy,  and 
that  the  presence  of  the  resin  reach  layer  did  not  significantly  affect  the  prediction. 

Evaluation  of  damage  distribution  shows  qualitative  agreement  with  the  radiographic  image. 
The  residual  stresses  were  not  accounted  for  in  the  analysis. 

An  internal  state  variable  approach  to  establish  property  degradation  rules  was  taken 
in  [24,25],  Two  failure  mechanisms  were  incorporated  in  [24]:  the  matrix  transverse 
cracking  and  fiber  tensile  breakage.  The  degraded  stiffness  matrix  due  to  transverse  cracking 
was  represented  as 

C  =  C?.eV/ 

ij  IJ 

where  C°  was  the  initial  (before  damage)  stiffness  matrix,  a/ the  damage  variable,  and  ky  the 

rates  of  individual  stiffness  component  degradation,  established  from  the  self-consistent 
model.  Fiber  breakage  was  detected  from  maximum  strain  failure  criterion.  Commercial 
finite-element  package  SAMSEF  was  chosen  to  implement  the  methodology.  Numerical 
results  presented  were  for  the  same  laminates  as  in  two-dimensional  analysis  [17],  It  was 
concluded  that  more  realistic  transverse  cracking  patterns  in  the  off-axis  plies  were  observed. 
The  loading  was  modeled  only  up  to  about  85  percent  of  the  failure,  due  to  numerical 
concerns.  No  residual  stress  or  delamination  effects  were  considered. 

The  discrete  internal  state  variable  approach  developed  in  [17]  for  two-dimensional 
analysis  was  extended  in  [25]  for  three-dimensional  prediction  of  bearing  strength  in 
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composite  bolted  joints.  A  rigid  bolt  was  modeled  and  ABAQUS  finite  element  package  was 
used  to  facilitate  the  property  degradation  model.  Hashin’s  failure  criterion  [12]  was  used  to 
identify  the  damage  mode  in  the  element.  Two  different  sets  of  internal  state  variables,  Df 
and  Df,  associated  with  tensile  and  compressive  loads  were  used.  Linear  eight-node 
elements  with  reduced  integration  were  utilized  for  the  three-dimensional  modeling  with  one 
element  through  the  thickness  of  the  ply.  The  delaminations  and  residual  stresses  were 
neglected.  In  the  bearing  failure  mode,  the  analysis  was  stopped  after  the  damage  propagated 
outside  the  washer-covered  area  around  the  hole.  The  model  proved  capable  of  predicting 
the  mode  of  failure:  tensile,  shear-out  or  bearing,  and  the  initial  load  drop-off.  No  detailed 
correlation  of  the  damage  progression  prediction  was  made  with  the  experiment. 

A  complex  boundary  condition  contact  problem  of  transverse  bending  of  a  composite 
plate  with  a  countersunk  bolt  was  considered  in  [26],  Complete  element  failure  approach 
was  adopted  for  property  degradation  modeling.  In  this  model  the  element  stiffness  was 
completely  degraded  if  any  one  of  the  incorporated  failure  criteria  was  met  for  the  average 
stress  component  values  of  the  element.  In  [26]  the  maximum  stress  failure  criterion  was 
used  for  in-plane  damage  and  Ye  [27]  delamination  failure  criterion  for  damage  caused  by 
out-of-plane  stresses.  Thus  through  the  complete  element  failure  approach,  the 
delaminations  were  included  in  the  model.  No  residual  stresses  were  considered.  The  model 
was  shown  to  produce  load  deflection  curves  in  good  agreement  with  the  experimental 
results  for  the  [90g/0s]s  and  [08/90x]s  laminates.  However,  no  detailed  correlation  of  the 
damage  progression  prediction  was  made  with  the  experiment. 

Both  two-  and  three-dimensional  models  have  been  developed  for  strength  prediction 
of  composites  with  open  and  filled  holes.  The  property  degradation  methodologies  used  for 
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three-dimensional  analysis  were  direct  extensions  of  the  ones  developed  for  two-dimensional 
analysis.  The  multi-scale  nature  of  the  failure  process  in  composites  requires  certain 
assumptions  to  be  made  in  order  to  formulate  the  property  degradation  algorithm.  The 
empirical  constants  entering  the  analysis  provide  tuning  mechanisms  in  order  to  achieve 
accurate  strength  prediction  for  a  given  material  system  within  a  certain  range  of  parameters. 
This  was  demonstrated  in  numerous  two-dimensional  studies  reviewed  above.  Note  that  only 
one  method  [17]  takes  into  account  the  residual  stresses  in  the  laminates.  Three-dimensional 
progressive  damage  analysis  alleviates  one  level,  namely,  the  through-the-thickness  of  the 
laminate,  level  of  averaging  in  the  stress  analysis.  The  price  for  that  is  at  least  an  order  of 
magnitude  increase  in  computer  time  and  resources  for  an  entry-level  practical  laminate,  e.g. 
an  eight-ply  quasi-isotropic  symmetric  laminate.  The  redemption  for  this  treatment  should 
be  no  less  than  a  physically  meaningful  and  verified  damage  origination  and  progression 
sequence  prediction,  in  different  types  of  laminates  such  as  soft  and  0°  dominated  laminates 
made  of  material  systems  with  different  toughness. 

In  the  present  report  a  three-dimensional  property  degradation  algorithm  will  be 
implemented  to  account  for  the  residual  stress  and  delamination  effects.  The  focus  of  the 
study  will  be  on  investigating  the  mesh  type  and  density  effect  aspects  of  modeling  which  are 
specifically  important  in  three-dimensional  property  degradation-based  damage  initiation  and 
progression  modeling. 

2.1  Problem  Formulation 

Consider  a  rectangular  TV-layer  laminate  built  of  orthotropic  layers  with  length  L  in 
the  x-direction,  width  A  in  the  y-direction,  and  thickness  H.  Individual  ply  thicknesses  are 
hp=z(p,-zl|>1  where  z=2*p>  and  z=2*t>1)  are  upper  and  lower  surfaces  of  the  p-th  ply, 
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respectively.  The  origin  of  the  x,y,z  coordinate  system  is  in  the  lower  left  comer  of  the  plate, 
as  shown  in  Figure  1.  A  circular  opening  of  diameter  D  with  the  center  at  x=xc  and  y=yc  is 
considered.  Uniaxial  loading  is  applied  via  displacement  boundary  conditions  at  the  lateral 
sides  (x=0,L),  while  all  other  boundaries  were  kept  free: 


-  ux(0,y,z)  =  ux(L,y,z)  =  uL  / 2, 
Uy  (0,  y,z)  =  Uy  (L,  y,  z)  =  0, 


(11) 


The  transverse  displacement  is  not  constrained  aside  from  a  rigid  body  constant.  The 
constitutive  relations  in  each  ply  are  as  follows: 


^  =C^(£U  -ocffAT) 


(12) 


where  Cfjkl  and  ajj  are  elastic  moduli  and  thermal  expansion  coefficients  of  the  p-th 
orthotropic  ply,  and  AT  is  the  temperature  change. 


2.2  Variational  Considerations 

Consider  volume  V  and  sub-volume  Vi,  which  for  correctness  may  be  considered 
surrounding  a  crack  (Figure  6).  We  shall  impose  a  number  of  stress  conditions  to  be  satisfied 
within  Vj : 

Fmioij)  =  0,  m  =  1,2,...  .  (13) 

The  Reissner-Hellinger  variational  principle 

»  =  111 Kfe,  +  «;,-)>  +  U(e„))dV  -  jj  T,u,ds , 

r  1  s, 

where  T,  are  applied  external  tractions,  will  be  considered  here  and  below  summation  upon 
repeated  indexes  is  implied.  The  Lagrangian  multiplier  method  is  used  to  impose  additional 
conditions  (13): 
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Figure  6.  Additional  Stress  Conditions  Applied  in  the  Subvolume 
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The  independent  variations  of  (14)  yield: 


8u  :  <7;-  ■  =  0,xeV,  a ■  ■n  .  -T.,xe 
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where 


1  , 

%/>  =2(,/ 

As  follows  from  Equation  (15),  the  additional  stress  conditions  (13)  are  contributing  to 
displacement  strain  relationship  in  the  sub-volume  Vi.  For  some  practically  meaningful 
cases,  the  Lagrangian  multipliers  can  be  explicitly  obtained.  Let 
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Fn,=°qm’m  =  W 

For  q- 2  it  corresponds  to  a  crack  in  the  plane  X2=0,  as  shown  in  Figure  6.  Substituting 
Equation  (15)  into  the  second  equation  (16),  one  will  obtain  the  following  relationship  within 
volume  Vi 

&ij  =  C  ijkl(U(k,I)  ~~  ekl  ~  ekl  (17) 

where  8  k  is  Kronechker’s  delta.  First  assume  that  only  one  component  of  stress  <j  =0, 
then,  according  to  Equation  ( 1 7) 

^  _  C qpkAU(k,l)  ~  eki  ) 

~~  c  ’ 

1P1P 

where  summation  is  carried  out  only  upon  k  and  /.  Substituting  this  equation  back  in  (17) 
one  obtains 

®ij  =  C ijkl(U(k,l )  ~~  ekl  )  >  (18) 

and 

C  C 

r*  ~r  wp  ipki  /io\ 

~  L:;U  C  '  ^ ' 

^  qpqp 

It  follows  from  equation  (19)  that  the  degraded  stiffness  matrix  C*jk,  is  a  symmetric  matrix 

with  zeroes  in  the  row  and  column  qp.  In  order  to  satisfy  all  three  conditions,  <Jqj  =  0, 

j=l,2,3,  the  degradation  operator  (19)  can  be  applied  recursively.  It  can  be  shown  that  the 
final  degraded  stiffness  matrix  is  independent  of  the  order  in  which  (19)  is  applied.  Similar 
results  can  be  obtained  intuitively  by  degrading  Young’s  moduli  E2  and  G12  with  additional 
and  less  intuitive  considerations  regarding  Poisson’s  ratio  coefficients. 
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2.3  Progressive  Damage  Modeling 

In  all  further  analyses  the  ply  material  axes  are  oriented  as  follows:  xi  is  the  fiber 
direction,  and  the  X3  axis  is  directed  perpendicular  to  the  lamina  plane.  Three  types  of 
damage  will  be  implemented  in  the  analysis:  (i)  fiber  breakage  (<r, .  =  0 ,  j=l,2,3),  (ii)  in¬ 
plane  transverse  cracking  (<t2  .  =  0,  j— 1,2,3),  and  (iii)  interlaminar  transverse  cracking 
(<r3 .  =  0 ,  j= 1,2,3).  A  combination  of  these  damage  types  is  modeled  by  requiring  the 

combined  set  of  stress  components  to  be  zero.  Maximum  stress  failure  criterion  was  utilized 
to  sample  each  element  at  the  center.  Nine  coefficients  were  computed: 


1. 

2. 

3. 

4. 

5. 

6. 

7. 

8. 

9. 


Coefficients  f  are  always  positive  and f  =1  indicates  local  failure  mode.  Criteria  1  and  2 
indicate  failure  breakage  and  lead  to  conditions  (i);  criteria  3,  4  and  9  indicate  in-plane 
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transverse  cracking  and  lead  to  conditions  (ii);  criteria  4-8  indicate  interlaminar  failure  and 
are  modeled  by  applying  stress  conditions  (iii). 

Progressive  damage  modeling  is  performed  by  incrementing  the  damage  rather  then 
the  load.  Each  step  of  this  procedure  calculates  the  load  increment  required  to  introduce 
additional  damage.  In  case  of  linear  constitutive  law,  the  first  load  increment  obtained 
through  this  scheme  corresponds  to  damage  initiation.  To  take  into  account  the  residual 
stress,  two  boundary  condition  problems  with  different  boundary  conditions  are  solved  at 
each  step.  One  of  them  is  the  mechanical  uniaxial  loading  (1 1)  under  unit  load  with  zero 
thermal  stresses,  and  the  other  is  the  thermal  loading  with  all  free-edge  (less  rigid  body 
motion)  boundary  conditions,  corresponding  to  residual  stresses.  The  flow-chart  of  the 
progressive  damage  analysis  algorithm  is  as  follows: 


1 .  Input  and  zero  damaged  elements. 

2.  Solve  the  mechanical  loading  problem  (11)  at  ul=0.001L  for  AT=0  and  determine  //w  in 
each  element. 

3.  Solve  the  thermal  loading  problem  for  given  AT  under  all  free-edge  boundary  conditions 
and  determine  //  in  each  element. 


4. 


Calculate  loading  coefficient 


P  = 


min 

j 

all  elements 


and  record  the  element  and  mode  for 


which  it  was  calculated. 

5.  Determine  damage  mode  (I),  (ii),  (iii)  from  the  previous  step.  Degrade  the  properties  by 
consecutively  applying  degradation  operator  (19). 

6.  Go  to  step  2. 
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The  procedure  is  implemented  in  computer  code  SVELT1 .2  allowing  for  various  types  of 
displacement  approximation  basis  functions. 

2.4  Numerical  Results 

Numerical  results  and  experimental  data  were  obtained  for  a  material  system 
IM7/5250-4  by  using  the  set  of  elastic  properties  described  in  Table  3.  In  all  cases  linear 
displacement  approximation,  corresponding  to  eight-node  brick  elements,  was  used,  with  the 
exception  of  the  geometry  models.  In  the  present  work  exact  curvilinear  transformations 
were  used  as  opposed  to  isoparametric  elements.  The  strength  properties  in  fiber  and 
transverse  directions  were  obtained  by  testing  unidirectional  eight-ply  thick  laminates  and  the 
shear  properties  by  testing  [45/-45]4S  laminates.  All  specimens  were  7”  wide  and  2”  long. 

The  notched  specimens  had  either  a  0.25-inch  or  0.5-inch  central  hole.  Unnotched  quasi¬ 
isotropic  [45/0/-45/90]  s  coupons  were  analyzed  first.  The  predicted  and  experimental  stress- 
strain  curves  are  shown  in  Figure  7.  The  beginning  of  the  numerical  curve  corresponds  to 
failure  initiation  and  correlates  well  with  experimental  values.  The  analysis  predicts  that  the 
failure  initiation  load  is  very  close  to  the  90°  ply  failure  load  [all  elements  fail  in  mode  (ii)]. 
At  approximately  0.77  percent  strain,  a  stiffness  loss  is  predicted  which  is  due  to  transverse 
cracking  of  the  ±45  plies.  The  ultimate  strain  is  controlled  by  fiber  failure  and  is  the  same 
for  experiment  and  prediction.  The  significant  strength  undeiprediction  is  essentially  due  to 
overprediction  of  the  stiffness  degradation  in  the  ±45  plies.  The  over-softening  effect  of  the 
present  analysis  is  due  to  complete  stiffness  degradation  in  the  transverse  direction  due  to 
transverse  cracking  (ii).  The  remedy  in  this  case  will  be  limiting  the  transverse  stiffness 
degradation  by  lower  limit  values  due  to  saturation  of  the  transverse  cracking  density 
[29,30]. 
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Axial  stress  [psi] 


Table  3 

Material  and  Strength  Properties  Used  in  Modeling 


1  material  type  (1 -elastic  orthotr.)  [psi] 


El  1=2.2 15E+07 

E22=1.300E+06 

E33=1.300E+06 

U13=3.270E-01 

U23=4.770E-01 

U12=3.270E-01 

G13=8.670E+05 

G23=4.770E+05 

G12=8.670E+05 

A1  l=2.500E-07 

A22=1.400E-05 


Xt=3.497E+05 

Yt=9.610E+03 


S=1.451E+04 


Xc=3.031E+05 

Yc=4.139E+04 


S 1 3=1 .25 1E+04 


Stress  free  temperature  T=350°F 


sx  [45/0/-45/90]s  -test 


Axial  strain,  % 

Figure  7.  Experimental  and  Predicted  Stress-Strain  Curves  of  a  [45/0/-45/90]s 
Laminate. 
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Three-dimensional  damage  propagation  modeling  in  notched  laminates  imposes  severe 
requirements  on  the  capability  of  the  finite-element  modeling.  Namely,  it  is  natural  to  expect 
that  a  three-dimensional  finite-element  model,  which  is  successful  in  predicting  the  strength 
of  the  laminate,  should  be  able  to  predict  the  strength  of  each  of  its  plies.  Thus  the  strength 
prediction  in  the  [Ox  ]  laminate  is  considered  next.  The  experimentally  observed  failure  mode 
of  this  laminate  is  development  of  transverse  cracks  or  splits  emanating  from  the  hole  which 
reduce  the  stress  concentration  in  the  fiber  direction,  and  provide  the  strength  of  the  laminate 
equal  to  that  of  the  ligament  strength.  This  mechanism  will  be  present  in  the  laminate,  where 
the  splits  may  not  propagate  through  the  entire  length  of  the  laminate  and  will  only  provide 
partial  stress  concentration  relief.  Typical  finite-element  subdivision,  shown  in  Figure  8,  will 
be  utilized.  The  progressive  failure  analysis,  however,  predicts  early  onset  of  fiber  breakage 
due  to  inability  to  model  the  stress  relief  in  the  fiber  direction.  To  illustrate  this  fact  the 
following  numerical  examples  were  considered.  In  the  finite-element  mesh  shown  in  Figures 
8a  and  8b,  all  the  elements  with  the  center  coordinates  Xj,yi,z,- ,  such  that  yc  -  D/2<  y,  <  yc  + 
D/2,  were  assigned  to  have  in-plane  matrix  cracking-type  damage  (ii).  This  band  of  elements 
is  darkened  in  Figure  8.  The  mesh  size  is  chosen  sufficiently  small,  so  that  the  roughness  of 
the  edges  of  this  band  near  the  hole  is  much  smaller  than  the  hole  diameter.  The  radial 
distribution  of  the  hoop  stress  (Txt  /  <70  normalized  to  applied  far-field  stress  is  shown  in 
Figure  9.  The  two  mesh  densities  yield  practically  equal  hoop  distribution  for  virgin 
laminates.  In  the  presence  of  massive  transverse  damage,  however,  the  hoop  stress  is 
reduced  very  little.  Significant  relief  of  the  stress  concentration  can  only  be  observed  if  the 
stiffness  of  the  selected  finite  elements  is  completely  reduced,  i.e.  (i)+(ii)  degradation 
conditions  applied.  This  assumption  is  physically  incorrect  since  no  massive  fiber  breakage 
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(b) 

Figure  8.  Radial  Type  Mesh  Configurations  Used  for  Predicting  the  Fiber  Stress 

Relaxation  due  to  Transverse  Splitting:  (a)  Coarse  Mesh;  (b)  Dense  Mesh. 
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Figure  9.  Normalized  Hoop  Stress  Oee/co  Distribution  as  a  Function  of  the  Distance 
from  the  Hole  Edge  in  the  Direction  Perpendicular  to  Loading. 
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in  the  selected  regions  is  observed.  The  effect,  pointed  out  above,  appears  to  be  finite- 
element  mesh  type  dependent.  A  different  type  of  mesh,  with  the  element  edges  aligned  with 
the  damage  propagation  direction,  was  used  next  for  analysis.  The  mesh  is  shown  in  Figure 
10.  Similar  tests  conducted  with  this  mesh  showed  that  by  applying  the  transverse  cracking- 
type  property  degradation  to  the  elements  within  the  selected  band,  one  observes  a  complete 
relaxation  of  the  hoop  stress,  as  shown  in  Figure  9.  Progressive  damage  modeling  performed 
with  this  mesh  for  unidirectional  laminates  with  two  hole  sizes  is  shown  in  Figure  11.  The 
experimentally-obtained  strength  value  for  the  0.25-inch-diameter  hole  is  307  ksi  and  is  very 
close  to  the  prediction. 

The  new  mesh  type  does  not  represent  a  solution  to  the  mesh  type  dependent  behavior 
of  damage  propagation  and  stress  redistribution.  There  are  no  reasons  to  expect  that  the 
mesh  shown  in  Figure  10  will  be  more  suitable  for  progressive  damage  analysis  in  off-axis 
plies  than  the  radial  mesh.  It  illustrates,  however,  the  additional  difficulties  in  three- 
dimensional  property  degradation  modeling  which  were  not  addressed  in  previous  works 
[23-26]. 

2.5  Conclusions 

A  3-D  progressive  damage  analysis  method  has  been  developed.  Nine  maximum 
stress  failure  criteria  have  been  used  for  local  failure  prediction.  Residual  stresses  have  been 
included  into  the  analysis. 

Damage  propagation  modeling  in  notched  specimens  shows  significant  dependence 
upon  mesh  type  and  density. 

Accurate  transverse  cracking  prediction  (including  thickness  effect)  has  a  significant 
influence  on  final  failure  prediction  of  both  notched  and  unnotched  specimens. 
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Figure  10.  Configuration  of  the  Mesh  with  the  Element  Edges  Aligned  with  the 
Damage  Propagation  Directions. 
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3.  DAMAGE  INITIATION  AND  PROGRESSION  IN  MULTIDIRECTIONAL 

LAMINATES  WITH  A  HOLE 


3.1  Experiment 

A  series  of  tension  tests  were  conducted  to  document  the  initiation  and  propagation  of 
damage  in  multidirectional  composite  laminates  with  a  centrally  located  circular  hole. 

Testing  was  also  conducted  for  unnotched  tensile  coupons  to  determine  the  first-ply  failure 
and  ultimate  strength  to  compare  with  the  hole  specimens.  The  hole  specimens  were  tested 
under  a  series  of  incremental  loadings  to  document  the  initiation  and  progression  of  damage 
and  the  corresponding  strain  response.  Acoustic  emission  data  were  also  collected  during 
loading  and  compared  to  strain  measurements  and  radiographic  images  to  identify  significant 
damage  events.  Additional  analysis  was  conducted  by  comparing  experimental  strain  data  to 
results  obtained  from  elasticity  theory  and  from  spline  variational  theory. 

The  material  system  chosen  for  this  work  is  graphite/toughened  bismaleimide, 
IM7/5250-4.  Two  types  of  laminates,  [02/902]s  and  [±30/902]s,  were  considered.  The 
thermoelastic  properties  necessary  for  analysis  were  determined  from  testing  of  [(%r,  [90]  si 
and  [±45]2s-  Each  laminate  was  cured  according  to  the  manufacturer’s  recommended  cure 
cycle  and  was  postcured  at  440°F  for  a  duration  of  five  hours  and  40  minutes.  The  average 
fiber  volume  was  found  to  be  62  percent,  and  the  void  content  was  determined  to  be  less  than 
one  percent. 

In  multidirectional  laminates,  the  holes  were  drilled  using  a  carbide-tipped  steel  drill 
bit  with  an  operating  speed  of  800  rpm  based  on  an  optimization  study  to  minimize  drilling- 
induced  damage.  Prior  to  initial  loading,  a  radiographic  inspection  of  each  specimen 
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